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REVIEW 

Using Network Biology to Bridge Pharmacokinetics and 
Pharmacodynamics in Oncology 

DC Kirouac^ and MD Onsum^ 

If mathematical modeling is to be used effectively in cancer drug development, future models must take into account both the 
mechanistic details of cellular signal transduction networks and the pharmacokinetics (PK) of drugs used to inhibit their onco- 
genic activity. In this perspective, we present an approach to building multiscale models that capture systems-level architectural 
features of oncogenic signaling networks, and describe how these models can be used to design combination therapies and 
identify predictive biomarkers in silico. 
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THE CHALLENGE 

The last 60 years have witnessed an explosion in scien- 
tific data about the molecular underpinnings of health and 
disease. Unfortunately, this increase in information has not 
translated into an increase in pharmaceutical research pro- 
ductivity-over the same 60-year period the annual output 
of new drugs has held roughly steady while R&D expenses 
have grown exponentially.^ Why is this? Some have specu- 
lated that all of the "low-hanging therapeutic fruit" has been 
harvested, and others have opined that the increasingly 
competitive market and regulatory environment is demand- 
ing more stringent benefit/risk profiles for new therapies.^ 
While these issues may hold some truth, what has become 
increasingly clear is that the development of future medicines 
will require a deeper appreciation of the complexity of dis- 
ease states, and an understanding of how drug treatments 
affect not only diseased tissue, but also surrounding healthy 
tissue and organs across diverse patient populations. One 
proposed solution to improve drug R&D productivity is to uti- 
lize computational modeling and simulation as a foundational 
platform, as it has been in other high-tech industries. By pro- 
viding a coherent framework for integrating knowledge and 
formalizing assumptions, mathematical models can enhance 
our ability to make informed and quantitative decisions. This 
is not a new idea - the FDA advocated "model-based drug 
development" as a pillar of its critical path initiative in 2004.^ 
The question is: why are notable successes on this front so 
scarce a decade later? 

A simplistic explanation is that biology is too complex- 
more complex than we comprehend, and still too complex 
to be formulated mathematically. While we may lack a uni- 
fied theory of molecular biology, this doesn't necessarily 
preclude our ability to develop models capable of answering 
specific questions, within a defined scope." Computational 
modeling has in fact played an important role in biomedical 
R&D for years. The issue is that the types of models cur- 
rently employed exist in parallel universes, with different 
academic journals, university departments, and languages.'' 



Computational systems biology is largely an academic disci- 
pline, focused on discovery-stage research with an empha- 
sis on understanding molecular and cellular mechanisms. 
Given the complexity of molecular biochemistry in compari- 
son with the sparse data typically available, such models are 
by nature generally nonidentifiable.'' At the other end of the 
spectrum, pharmacometrics is focused on utilizing available 
clinical data to guide pressing drug development decisions. 
The technical emphasis is thus on robust parameterization 
and predictive accuracy rather than mechanistic detail, which 
biases towards more empirical model formulations. Develop- 
ing multiscale models both with respect to the underlying sci- 
ence (connecting genes and proteins to cellular phenotypes, 
to tissue, organ, and population-level responses) and the 
drug development process (linking discovery through clinical 
development) is a critical challenge.^ While a number of orga- 
nizational and cultural issues are involved, a key technical 
difficulty lies in linking together the different time and length 
scales at play into a single representative format. 

SIGNALING NETWORKS LINK MOLECULAR EVENTS 
TO CELLULAR RESPONSES 

Signal transduction networks function at the interface of the 
cell's external and internal worlds, converting environmental 
information into biochemical changes, and ultimately discrete 
cell fate decisions. These molecular circuits thus serve not 
just as linear messengers, but as computational systems 
performing highly complex, multivariate signal processing. 
Many pathologies, in particular cancers, are currently under- 
stood to arise from aberrant cellular information processing, 
a result of genetic and epigenetic perturbations in signaling 
networks.'' Signal transduction networks thus appear to be 
an appropriate and perhaps necessary means to mecha- 
nistically link PK and pharmacodynamic (PD) models." Tra- 
ditional drug exposure-response models are constructed 
using empirical transfer functions, though often including 
semi-mechanistic details based on receptor theory.^" These 
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are generally sufficient for purposes in which they were 
developed — defining dosing regimens for monotherapies, 
based on existing clinical data. Such empirical formulations 
however fall short of the most pressing challenges in modern 
drug development. Using oncology as an example, where an 
increasing array of molecular targeted therapies are available 
for clinical use, meaningful responses to single agents are 
rare and generally transient, necessitating combination drug 
regimens." In addition, tumors are widely heterogeneous at 
the molecular level. This molecular diversity mediates differ- 
ential drug sensitivities between patients, and even between 
clonal cell populations within individual tumors. Predic- 
tive molecular biomarkers are thus required to match drug 
regimens with patient populations expected to respond to 
therapy. Testing drug combination regimens in an empirical 
manner is, however, becoming increasingly unfeasible. Given 
the number of anticancer agents currently available and the 
genetic diversity of the disease, the combinatorial explosion 
of possibilities is beyond what can typically be tested through 
brute force screening. Computational modeling and simula- 
tion could provide means to beset this limitation, allowing 
one to perform ultra-high throughput in silico screening and 
analysis of synthetic treatment-response data. Developing 
mathematical models capable of meeting this challenge is, 
however, a nontrivial task. Below we describe the obstacles 
and outline some proposed solutions. 

QUANTITATIVE LOGIC AS A FRAMEWORK FOR 
REPRESENTING SIGNALING NETWORKS 

So how does one mathematically formulate network biology- 
based PK-PD models? Most biochemical network models 
published to date have been based on mass action kinetic- 
based ordinary differential equations. Motivation for doing 
so arises from the desire to develop mathematical models 
based on fundamental physiochemical properties and reac- 
tion constants, thus transferable across different cells, tis- 
sues, or disease states. Moreover, given the extensive history 
of mass action kinetic-based modeling in other disciplines, 
established methods and expertise are widely available to 
draw from. There are many notable successes, both in funda- 
mental cell biology"'^^ as well industry pursuits such as drug 
target discovery"' and therapeutic antibody design. ^^-^^ How- 
ever, caution should be used when applying the assumptions 
underlying mass action kinetics to intracellular processes. 
Most biochemical reactions involved in cellular signal trans- 
duction take place as part of multiprotein complexes, often 
tethered to scaffolds or cell membranes. The kinetics would 
therefore be expected to deviate from that predicted by the 
laws of mass action, which assume homogeneous solution- 
phase reactions. More importantly, as our knowledge of 
molecular biology is still far from complete, it is very difficult 
to parse biochemical cascades down to fundamental reaction 
steps, or to account for all relevant molecular species and 
reactions.^" This problem is particularly acute for processes 
downstream of canonical signaling cascades, connecting 
signaling events to gene expression changes or cellular phe- 
notypes. Even extensively detailed physiochemical-based 
models thus often contain many lumped parameters, which 
must be estimated by fitting to experimental data rather than 



derived from biophysical properties. This may underlie one 
of the difficulties in extrapolating model parameters across 
different cell lines. 

It is also important to recognize the distinct time scales 
at play. Dynamic events triggered by cell surface recep- 
tor engagement reach (quasi)-steady state within minutes 
to a few hours, while phenotypic readouts (i.e., measurable 
changes in bulk tumor size) are typically quantified on the 
order of days to weeks. As a consequence, many molecu- 
lar events may be represented algebraically rather than with 
more arduous differential equations. Another practical con- 
sideration is the type of data available for model training. 
Biochemical measurements are typically semi-quantitative, 
lacking the precision (molecules/cell) and the coverage (mea- 
sured vs. inferred species) required to parameterize mass- 
action kinetic-based models. Quantitative logic provides an 
alternative and relatively simple formalism to represent the 
structure and information processing capabilities of signaling 
networks, ^^'^^ bridging the distinct time scales of biochemi- 
cal and physiological events. Quantitative logic networks 
are assembled using Hill-type equations, malleable signal- 
response curves representing information flow between 
nodes (i.e., protein species). When a network node contains 
multiple inputs, quantitative logic gates can be used to repre- 
sent various types of signal processing. These are analogue 
extensions of Boolean logic truth tables, the most common 
forms being AND, NAND, OR, and NOR gates, which can be 
configured to recapitulate biochemical and pharmacological 
mechanisms. The algebraic equations can be easily extended 
into differential equation form, so as to capture both fast 
(steady state) and slower (dynamic) process together, using 
systems of differential-algebraic equations (Figure 1). The 
logic gates and hill functions used to describe signal flow in 
quantitative logic networks are data-driven rather than based 
upon fundamental biophysical constants. However, they are 
in fact not as different from mass action kinetic ordinary dif- 
ferential equations as initial appearances suggest, given that 
such models often contain many data-driven parameters as 
well. Note there is no single best approach to modeling cell 
signal transduction. The choice between alternatives, from 
purely data-driven statistical models to physiochemical ODEs, 
should be determined by the specific questions at hand, data 
available, and specifics of the underlying biology. 

The next challenge is specifying the level of detail required 
to describe the relevant biology and mechanistically capture 
exposure-response relationships. Explicitly characterizing 
all biochemical reactions taking place in a cell ("bottom-up" 
network reconstruction) is currently impractical given the 
complexity and remaining knowledge gaps involved in cel- 
lular signal transduction.^" Can we then identify archetypal 
properties of signaling networks necessary to be included 
in our models, striking a balance between predictive power 
and mechanistic representation of the biology? Below we 
describe four systems-level design features as hallmarks of 
oncogenic signaling networks: modularity, redundancy, adap- 
tation, and heterogeneity. We will outline how these features 
can be captured mathematically, and the implications for 
drug development. While our focus is oncology, we believe 
the principles and approaches described are generalizable 
to other areas. 
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Figure 1 Mathematical underpinnings of quantitative logic networks, (a) Hill-type equations can be used to describe input-output relationships 
between individual nodes (i.e., proteins) in a network. Using two parameters, the Hill coefficient (K) and EC50 value (t), these malleable 
equations are capable of describing a wide range of stimulatory and inhibitory response curves, (b) Hill-type equations can be modified to 
capture dynamic as well as algebraic (steady state) relationships. Shown are exemplary time course phase portraits of an output variable (Y) 
in response to a pulse of stimulator (X), indicated by dashed lines. The system gain (G.^) and symmetric rates of production/degradation (r^) 
can be tuned to capture differing system dynamics, while the basal set-point of Y (Yg) quantifies the steady-state level, (c) Various types of 
quantitative logic gates can be used to capture biologically relevant responses to dual inputs. Shown are four exemplary surface responses: 
dual stimulatory OR gate (i.e., alternative RTKs activating PI3K, with different relative activation weights (w)), dual stimulatory AND gate (i.e., 
co-dependent transcription factors activating gene expression), competitive inhibition (i.e., stimulatory ligand vs. ligand-blocking antibody), 
and non-competitive inhibition (i.e., a stimulatory ligand vs. small molecule RTK inhibitor). PI3K, phosphatidylinositide 3-kinase. 



HALLMARK 1: MODULARITY 

It has been long standing practice in molecular biology to 
classify proteins and their biochemical interactions into 



pathways. It is now appreciated that such simple classification 
schemes, and original depictions of pathways as linear, insu- 
lated signaling cascades are vast oversimplifications of more 
complex integrated networks. However, evidence supports 
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the notion of pathways as functional networl< modules — clus- 
ters of nodes (proteins) which are highly interconnected but 
less connected to other parts of the network, and mapped 
to highly overlapping phenotypes.^" This concept is per- 
haps best supported by data emerging from cancer genome 
sequencing projects. Even within clinically homogenous dis- 
eases, mutations between tumors are found to be incredibly 
disparate. However, when mutations are mapped onto pro- 
tein interaction networks, they fall onto a relatively conserved 
set of oncogenic modules. Most notable are the phospha- 
tidylinositide 3-kinase (PI3K)/AKT and mitogen-activated 
protein kinase/extracellular signal-regulated kinase (IVIAPK/ 
ERK) cascades, and the receptor tyrosine kinases (RTKs) 
which regulate their activities, the TP53 and cell-cycle con- 
trol modules, the apoptotic caspase machinery, as well as a 
handful of other more indication and tissue-specific signal- 
ing cascades. Moreover, mutations within these network 
modules are generally found to be mutually exclusive (the 
prototypic example being P/SKC/A-activating mutations vs. 
P7E/\/ deletions), indicating that what "matters" to the cell is 
the activation of modules rather than individual molecules. 

The implication is that mathematical representations of 
molecular networks can be condensed to more abstract, 
modular formats, while maintaining functional relevance. This 
can be achieved using rigorous bioinformatic^'' or graph theo- 
retical approaches, or heuristically (as is common practice) 
by selecting key sentinel nodes representative of modules 
activation states (i.e., phosphorylated AKT and ERK for the 
RISK and MARK cascades, respectively). This reduction 
vastly simplifies the computational requirements for param- 
eterizing sparsely sampled and underdetermined systems. 
However, if the module decomposition scheme is appropri- 
ate, functional inferences from such granular descriptions 
should be equivalent to more detailed molecular forms. 

HALLMARK 2: REDUNDANCY 

Redundancy is observed at multiple levels within biochemi- 
cal networks - we have already mentioned one example in 
that mutations within a network modules display functionally 
equivalent (or at least highly similar) effects. This is a con- 
sequence of biochemical information flow - disrupting a sig- 
nal at different points along the transmission route will have 
similar effects ("vertical" redundancy). An alternative form of 
"horizontal" redundancy arises from cell's capacity for alter- 
native biochemical wiring. A consequence of evolutionary 
recycling of protein-binding domains, network modules and 
their constitutive proteins are capable of combining in alter- 
nate configurations. The family of RTKs for example, though 
expressed and employed in a functionally diverse manner, all 
transmit signals across the plasma membrane through bind- 
ing to a limited set of adaptor proteins. As a result, they inevi- 
tably plug into the same PI3K and MARK cascades, though 
with differing strengths and dynamics. The functional conse- 
quences of RTK activation are thus highly overlapping, as 
the signals converge on a limited set of core cytoplasmic hub 
proteins (so called "bow tie" architecture).™ Clinical evidence 
now supports this theory, as different RTKs have been found 
to be functionally exchangeable even within individual tumor 
biopsies.^" The consequences of this architectural feature are 



becoming increasingly appreciated in oncology through the 
phenomena of "pathway switching." Tumors displaying onco- 
gene addiction to a single RTK driver are often capable of 
switching dependency to an alternate, but functionally equiv- 
alent RTK following inhibitor treatment. 

The core cytosolic cascades themselves also display a 
degree of horizontal redundancy. For example, in a variety of 
preclinical cancer models it has been demonstrated that com- 
plete suppression of either PI3K or MARK pathways alone 
results in at best tumor growth suppression (or stasis), while 
cotargeting both cascades is required to induce apoptosis 
and tumor regression. This phenomenon can be described 
mathematically using logical OR gates along the signal trans- 
mission route from surface receptors to phenotypic readouts 
(i.e., AKT can be activated by EGFR or ERBB3 receptors, 
apoptosis is inhibited by active AKT ox ERK' signals). 

HALLMARK 3: ADAPTATION 

A distinguishing feature of living vs. nonliving matter is the 
ability to maintain internal stasis in the face of a continually 
shifting environment. Homeostatic regulation typically involves 
feedback control circuits, and many biochemical feedback cir- 
cuits have been characterized. Feedback controls within signal 
transduction cascades buffer the effect of targeted inhibitors 
and other perturbations, an example being phosphorylated 
ERK inhibiting the activation of its upstream kinases.^^ The 
activities of intracellular signaling cascades are also fine-tuned 
through regulation of cell surface receptor expression. ^''■^^ 
These receptor-coupled feedback circuits underlie the so called 
"Whac-a-Mole®" effect, named after the popular arcade game, 
wherein suppression of a single RTK or pathway induces com- 
pensatory activation of a parallel, functionally redundant path- 
way These feedback regulatory circuits can be simulated using 
methods adapted from engineering control theory, and inte- 
grated within logic-based signaling network models (Figure 2). 
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Figure 2 Adaptive feedback circuits based on control theory. Slight 
modification of the equation presented in Figure lb can be used 
to describe the activation of node (Y) by one of its downstream 
dependents (X), forming a negative feedback control circuit. Two 
alternative configurations are depicted, each displaying similar 
response dynamics to targeted inhibition (I) of X - up-regulation of 
Y to a new steady-state value. For exemplary purposes, the inhibitor 
(I) is constant in these simulations, but a pharmacokinetic driving 
model could be incorporated to account for more complex in vivo 
pharmacodynamics. 
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HALLMARK 4: HETEROGENEITY 

Whereas tumors are generally believed to be monoclonal 
in origin, full-blown carcinogenesis results from a brandied 
evolutionary process, resulting in extensive clonal heteroge- 
neity at diagnosis.^^'^"'^' The genetic and biochemical diver- 
sity of clones within a tumor then serves as the fodder of 
evolutionary adaptation. A result of the diversity in molecular 
wiring, biochemical responses to drug treatments would be 
expected to differ between clonal populations, underlying dif- 
ferential drug sensitivities.^'*'^'' Selective pressure then allows 
for the expansion of rare drug resistant clones and eventual 
clinical relapse."" This is perhaps the primary obstacle limit- 
ing our ability to achieve sustained remissions with targeted 
anticancer agents. For translational models to be clinically rel- 
evant, it is thus critical to understand molecular mechanisms 
of drug sensitivity and resistance, both between patients as 
well as between different clones within an individual tumor, 
and describe the phenomenon mathematically. 

So how do we then attempt this feat? A critical first step 
is to connect internal signaling states to tumor growth, or 
other relevant phenotypes. This is a nontrivial task given the 
multivariate, non-linear nature of cellular signal decoding. ''^■''^ 
However, by designing appropriate perturbation-response 
experiments, it is possible to define proliferation and apoptosis 
rates as functions of a few key intracellular effector molecules 
(i.e., phospho-AKT and -ERK). Next, given an understand- 
ing of the variability in key nodes (protein abundance) and 
parameters (biochemical rates), IVIonte Carlo simulations can 
be employed to simulate heterogeneity in tumor growth rates 
among molecular diverse tumor "clones." Performing these 
simulations in both the control and drug treatment conditions. 



one can then quantify the sensitivity of each clone to drug 
treatment using metrics such as percentage tumor growth 
inhibition. Such synthetic multivariate datasets can then be 
used to extract statistical relationships between biochemi- 
cal variability and drug responsiveness (Figure 3). Given the 
semi-empirical nature of quantitative logic networks, param- 
eter variations reflect relative changes from a prototypic cell 
type used for model training. Mutations can be represented 
by their functional effects. For example, gene amplification 
(or other means of increased protein expression) can be 
captured by simply increasing the relative basal level of a 
node, and conversely for gene knockouts. Activating kinase 
mutations can be captured by decreasing the EC50 value in 
upstream activating connections and thus simulating hyper- 
activity, conversely for deactivating mutations. While their 
parametric changes may be viewed as crude representations 
of molecular biochemistry, their relative effects on network- 
level responses to drug treatment may still be informative. 

Variations on pieces of this approach have been employed 
elsewhere. Population pharmacokinetic modeling is widely 
used to understand how interindividual variability in drug 
metabolism and physiology affects drug exposure and clini- 
cal responses."'' Tumor dynamic models have been used to 
predict overall survival based on the variability in initial tumor 
growth responses to drug treatment.'"'''"' At the cellular level, 
simulations of molecular variation have been used to recapit- 
ulate phenotypic heterogeneity in cell populations."" Signal- 
ing networks provide a means to link these different classes 
of models, connecting drug exposure, through molecular 
activities, to clinically relevant responses. This could enable 
the identification of mechanistic, rather than statistically- 
based clinical biomarkers for patient stratification. 
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Figure 3 Heterogeneity in tumor growth and drug responsiveness 
described using Monte Carlo simulations. Tumor growth models (in 
this case exponential, but other forms are possible) describe the bulk 
tumor burden overtime as a balance between population-average 
cell proliferation (|j) vs. death (5) rates. These phenotypic parameters 
can in turn be described as functions of cell signaling states. By 
randomly varying model parameters describing signal flow through 
the network (i.e., protein levels, biochemical kinetics, and gene 
mutations) and performing repeated Monte Carlo simulations, it 
is possible to create synthetically heterogeneous tumor "clones." 
Following in silico treatment of tumor populations with drugs or 
drug combinations, statistical analyses can extract mechanistic 
biomarkers predictive of drug sensitivity and resistance. 



STRATEGY AND FUTURE PERSPECTIVES 

Combined, these four characteristic features of oncogenic 
signaling networks underlie the robustness of tumors to tar- 
geted perturbations, limiting the efficacy of drug treatment. 
Cancer cells addicted to an oncogenic pathway often have 
multiple built-in escape routes leading to drug resistance. 
Quantifying the mechanisms by which signaling networks 
mediate drug sensitivity vs. resistance, however, could enable 
the design of synergistic drug combinations, and the identifi- 
cation of predictive stratification and PD biomarkers a priori. 
Note that the mathematical concepts we have described are 
borrowed from the scientific literature — the use of Boolean 
networks in molecular cell biology has an extensive history."^ 
What have been lacking until very recently were the experi- 
mental tools necessary to generate sufficiently quantitative 
and thorough data to parameterize these network models. 
The availability of highly selective kinase inhibitors against 
both cell surface receptors and intracellular transducers, cou- 
pled with quantitative, multiplex and high-throughput protein 
measurement technologies"" are finally providing the means 
to do so. Systematic perturbation-response experiments can 
be used to produce datasets necessary to specify the net- 
work connectivity, logic gates, and adaptive feedback circuits 
employed by cancer cells. By profiling phenotypic readouts 
(i.e., proliferation) in tandem with biochemical measure- 
ments, it is then possible to specify cell fate decisions as 
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Figure 4 Example representation of a generalized cancer cell 
signaling network. We have recently used a variation of this 
generalized network model for in silico combination drug screening 
in £flSS2-amplified breast cancer."" Three receptor tyrosine 
kinases (RTKs; in our specific case three ERBB-family receptors) 
are activated in ligand-dependent or independent manner. Active 
(phosphorylated forms) of the receptors initiate the canonical PI3K/ 
AKT and MAPK/ERK cascades, though with different characteristic 
biases, captured using multi-input OR gates. Active (phospho-) 
AKT and ERK in turn inhibit RTK expression via negative feedback 
regulatory circuits (red lines). Phospho-AKT and ERK were identified 
as driving cell proliferation and inhibiting cell death via logical OR 
and AND gates respectively. Measurable tumor cell growth then 
represents a balance between rates of cell proliferation and death. 
The signal transduction portion of the model (RTKs to AKT/ERK) 
was represented algebraically, while the RTK feedback circuits and 
cell-growth regulation were represented using differential equations. 
MAPK/ERK, mitogen-activated protein kinase/extracellular signal- 
regulated kinase; PI3K, phosphatidylinositide 3-kinase. 

multivariate functions of network activation states, a neces- 
sary step in developing translational models."" 

We have recently described such an approach to develop 
targeted combination therapies for the treatment of HER2- 
amplifed breast cancer,^" the general network schematic 
of which is summarized in Figure 4. A model connecting 
ErbB receptor family signaling to cell growth was param- 
eterized using in vitro profiling data, and in silico screening 
was then performed to identify synergistic drug combina- 
tions and predictive biomarkers. Select results were then 
validated in vivo using xenograft tumor models. ErbB3- 
mediated signaling was identified as a potent mechanism 
of resistance to multiple targeted agents, and combin- 
ing the ErbB3 inhibitor MIVI-IH with ErbB2 inhibitors was 
found to induce synergistic tumor regression. Other groups 
have recently published similar mechanistic PK-PD models 
(though parameterized entirely in vivo) for other preclinical 
cancer models and drugs. These include models of MEK 
inhibitor treatment in BRAF-mutant melanoma, rituximab 



treatment of non-Hodgkin's lymphoma, gefitinib treatment 
of glioblastoma.''^ 

Whereas this approach has proven fruitful in a few select 
cases, incorporating cellular networks into PK-PD models 
is by no means a panacea to drug development challenges. 
Computational models are only as good as their assumptions, 
and assumptions which hold in one context may break down 
in another. First, our knowledge about the topology of intra- 
cellular signaling networks remains far from complete. Even 
extensively studied signal transduction cascades, when recon- 
structed from unbiased proteomic methods turn have orders 
of magnitude more interactions than depicted in canonical 
signaling maps.^° Determining which of these canonical and 
noncanonical interactions to include in a model thus remains 
an open problem. More vexing, network connectivity appears 
to be fluid, such that biochemical circuitry differs between tis- 
sues, and even in the same cell under alternate environmental 
conditions. This poses a particularly acute challenge for trans- 
lating model predictions from cell culture, to animal models, 
to clinical strategies. Developing a clearer understanding of 
the variability in network architectures across different tissues 
and cancer indications will be essential to broadening the 
scope of mechanistic PK-PD models. Perhaps tumors can be 
clustered into distinct groups based on molecular wiring and 
pathway dependence, rather than histology or molecular pro- 
files? This would provide a more functionally relevant guide to 
drug development strategies. Another limitation of the above 
approach is that cancer cells are modeled as autonomous 
units, ignoring tissue architecture, physiology, and cell-cell 
interactions. ^"^^ These processes are known to play important 
roles in tumor biology and pharmacology, such that multiscale 
models must eventually describe tumors as complex tissues 
involving intercellular signaling networks,^** rather than inde- 
pendent collections of cells. 

Model development is an iterative process. Successive 
rounds of model parameterization, verification, and modifi- 
cation require a close collaboration between computational 
scientists, experimental biologists, and clinicians. Fostering 
such relationships can be challenging in large organizations 
with established cultures and clearly defined boundaries 
between departments and scientific disciplines. Overcoming 
such cultural barriers is perhaps just as challenging as the 
technical issues we have discussed. Despite the remaining 
obstacles, we believe the stepwise introduction of mecha- 
nistic, network biology-based mathematical models will be 
essential for translating knowledge about molecular biology 
into clinical decisions. The approach described and math- 
ematical tools presented represent our humble attempt to 
do so. Through the continued development and application 
of mechanistic model-based drug development, we hope to 
see a new generation of highly selective and effective drug 
regimens emerge to better serve the needs of patients. 
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